Ch 16 exercises

The aim of these exercises is to continue making sure you’re comfortable with handling multivariate data. In this chapter’s, you’ll focus on analyses based on dissimilarities/distances, including fitting linear models to these kinds of response variables.

For each of the examples below, you should follow the sequence we’ve used previously, as far as it’s sensible:

1. What is the biological question?

2. Is the predictor continuous or categorical?

3. Write out the linear model corresponding to this question.

4. What distribution do you expect the response variable to follow?

5. What are the assumptions behind the statistical model you’ll fit?

   1. Are those assumptions satisfied?

6. Fit the model

   1. How will you assess whether the model fits well?

   2. Can you detect an effect of the predictors?

   3. How do you measure the effect?

7. What do you conclude (including any cautions)

---

A

Dixon et al. (2018) focused on assemblages of reptiles in forests and woodlands of southeastern Australia, with a particular interest in relationships between these assemblages and the fire history of the landscape. They identified 81 sites that varied in time since fire, from 6 months to >96 y. Rather than a continuum, three categories of time since fire were used, 0.5-2y, 6-12y, and >96y. Sites were also classified according to habitat.

Reptile assemblages were sampled with a range of methods, including visual surveys and camera traps, to give counts of 20 reptiles, whose abundance ranged from 0-1 to 0-126, depending on species.

Data are available from dryad, as Rep_abund.csv. You’ll want to focus on the columns with reptile numbers, and the one with the fire history (tsf). For convenience, we’ve extracted those data as dixonbiota.csv

Start by having a quick look at the data.

df <- read_csv("data/dixonbiota.csv")

Rows: 79 Columns: 22── Column specification ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr  (2): site, tsf
dbl (20): adup, aplat, amur, amac, aram, dcor, ecun, esax, eul, htal, ldel, lgui, lwhi, ppor, pent, pspe, ptex, rdie, tnig, vros
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

head(df, 10)

The file is pretty straightforward, with tsf being the fire category, and columns 3-22 each representing one reptile taxon.

Outline how you’d assess whether the reptile assemblage (the combination of species present and their abundances) is related to fire history, using a dissimilarity-based approach.

You should think about how you’ll deal with the data:

• Will you need a transformation, standardization, etc.?

• What are the consequences of any decisions you make for the interpretation of your analysis?

• What measure(s) of dissimilarity are appropriate?

Run the analysis and provide your interpretation

MDS on reptile abundances

Start with dissimilarities

Students should have some discussion of raw vs transformed (e.g. 4th root), standardization.

Fourth root is common. Using it downweights the influence of abundant species, while raw allows those species to have a strong influence. There are questions about which one of these reflects ecological function better

Data are all abundances

df1 <- df[,-(1:2)]   #Stripped back file without labels, leaving only 
df1s <- wisconsin(df1)   #Common to do some standardisation
df1s.bc <- vegdist(df1s,'bray')   #With standardisation (wisconsin)

Visualize results

Need to decide whether to use 2 or 3 dimensions

#Run 3d first
df1s.mds3 <- metaMDS(df1s.bc,k=3,autotransform=FALSE,try=40,trymax=80,maxit=200)

Run 0 stress 0.184364 
Run 1 stress 0.1893708 
Run 2 stress 0.1873869 
Run 3 stress 0.1843639 
... New best solution
... Procrustes: rmse 0.00106849  max resid 0.004608248 
... Similar to previous best
Run 4 stress 0.185286 
Run 5 stress 0.183525 
... New best solution
... Procrustes: rmse 0.0237458  max resid 0.158325 
Run 6 stress 0.1843844 
Run 7 stress 0.1835265 
... Procrustes: rmse 0.001382192  max resid 0.01054027 
Run 8 stress 0.1853324 
Run 9 stress 0.1845136 
Run 10 stress 0.1835252 
... Procrustes: rmse 0.0001633524  max resid 0.0007231613 
... Similar to previous best
Run 11 stress 0.1838506 
... Procrustes: rmse 0.00923149  max resid 0.07077904 
Run 12 stress 0.1835255 
... Procrustes: rmse 0.00114584  max resid 0.008962794 
... Similar to previous best
Run 13 stress 0.1851463 
Run 14 stress 0.1835256 
... Procrustes: rmse 0.001019902  max resid 0.007445476 
... Similar to previous best
Run 15 stress 0.1835269 
... Procrustes: rmse 0.0006013782  max resid 0.004532644 
... Similar to previous best
Run 16 stress 0.1880208 
Run 17 stress 0.1835247 
... New best solution
... Procrustes: rmse 0.0002909292  max resid 0.002132903 
... Similar to previous best
Run 18 stress 0.1835254 
... Procrustes: rmse 0.0006796255  max resid 0.003265735 
... Similar to previous best
Run 19 stress 0.1880236 
Run 20 stress 0.1911568 
Run 21 stress 0.1893846 
Run 22 stress 0.1835257 
... Procrustes: rmse 0.0005354247  max resid 0.00385845 
... Similar to previous best
Run 23 stress 0.1901928 
Run 24 stress 0.1927536 
Run 25 stress 0.1894269 
Run 26 stress 0.1873883 
Run 27 stress 0.1844043 
Run 28 stress 0.1835262 
... Procrustes: rmse 0.0006881577  max resid 0.005441543 
... Similar to previous best
Run 29 stress 0.1959468 
Run 30 stress 0.1853514 
Run 31 stress 0.1888215 
Run 32 stress 0.1853335 
Run 33 stress 0.1835249 
... Procrustes: rmse 0.0006178072  max resid 0.003845607 
... Similar to previous best
Run 34 stress 0.1907431 
Run 35 stress 0.1859705 
Run 36 stress 0.187387 
Run 37 stress 0.1835254 
... Procrustes: rmse 0.0007605464  max resid 0.004125702 
... Similar to previous best
Run 38 stress 0.1907833 
Run 39 stress 0.1835253 
... Procrustes: rmse 0.0003930969  max resid 0.001570481 
... Similar to previous best
Run 40 stress 0.1875561 
*** Best solution repeated 7 times

stressplot(df1s.mds3, main="Shepard plot")

df1s.mds3

Call:
metaMDS(comm = df1s.bc, k = 3, try = 40, trymax = 80, autotransform = FALSE,      maxit = 200) 

global Multidimensional Scaling using monoMDS

Data:     df1s.bc 
Distance: bray 

Dimensions: 3 
Stress:     0.1835247 
Stress type 1, weak ties
Best solution was repeated 7 times in 40 tries
The best solution was from try 17 (random start)
Scaling: centring, PC rotation, halfchange scaling 
Species: scores missing

# Now look at 2d
df1s.mds2 <- metaMDS(df1s.bc,k=2,autotransform=FALSE,try=40,trymax=80,maxit=200)

Run 0 stress 0.2620287 
Run 1 stress 0.2653405 
Run 2 stress 0.2649768 
Run 3 stress 0.2609687 
... New best solution
... Procrustes: rmse 0.04264938  max resid 0.2810735 
Run 4 stress 0.2688938 
Run 5 stress 0.2677242 
Run 6 stress 0.2650354 
Run 7 stress 0.260286 
... New best solution
... Procrustes: rmse 0.03714916  max resid 0.2572555 
Run 8 stress 0.2628528 
Run 9 stress 0.2688503 
Run 10 stress 0.2760927 
Run 11 stress 0.261766 
Run 12 stress 0.2705515 
Run 13 stress 0.2641006 
Run 14 stress 0.2672022 
Run 15 stress 0.2624477 
Run 16 stress 0.2642981 
Run 17 stress 0.2760958 
Run 18 stress 0.2625464 
Run 19 stress 0.2616471 
Run 20 stress 0.2699231 
Run 21 stress 0.2557671 
... New best solution
... Procrustes: rmse 0.08490104  max resid 0.2915107 
Run 22 stress 0.292635 
Run 23 stress 0.2552758 
... New best solution
... Procrustes: rmse 0.01858635  max resid 0.105244 
Run 24 stress 0.2762298 
Run 25 stress 0.2712857 
Run 26 stress 0.2639198 
Run 27 stress 0.2705536 
Run 28 stress 0.2682175 
Run 29 stress 0.2638279 
Run 30 stress 0.2696993 
Run 31 stress 0.2650743 
Run 32 stress 0.2681181 
Run 33 stress 0.2647961 
Run 34 stress 0.2728345 
Run 35 stress 0.2705501 
Run 36 stress 0.2638431 
Run 37 stress 0.2602072 
Run 38 stress 0.2605345 
Run 39 stress 0.2700785 
Run 40 stress 0.2717778 
Run 41 stress 0.2657743 
Run 42 stress 0.2654575 
Run 43 stress 0.2689471 
Run 44 stress 0.2675718 
Run 45 stress 0.2591926 
Run 46 stress 0.2623515 
Run 47 stress 0.2624498 
Run 48 stress 0.2587334 
Run 49 stress 0.2738451 
Run 50 stress 0.2591982 
Run 51 stress 0.2643658 
Run 52 stress 0.2624203 
Run 53 stress 0.2666844 
Run 54 stress 0.2676336 
Run 55 stress 0.2665836 
Run 56 stress 0.2603051 
Run 57 stress 0.2624025 
Run 58 stress 0.2598046 
Run 59 stress 0.2567615 
Run 60 stress 0.2652927 
Run 61 stress 0.2670697 
Run 62 stress 0.2676958 
Run 63 stress 0.2652794 
Run 64 stress 0.2663321 
Run 65 stress 0.2624341 
Run 66 stress 0.2634596 
Run 67 stress 0.2761549 
Run 68 stress 0.261702 
Run 69 stress 0.2624381 
Run 70 stress 0.2618053 
Run 71 stress 0.2554227 
... Procrustes: rmse 0.01171501  max resid 0.0964735 
Run 72 stress 0.2655037 
Run 73 stress 0.2634403 
Run 74 stress 0.2820013 
Run 75 stress 0.2704852 
Run 76 stress 0.2620309 
Run 77 stress 0.2627632 
Run 78 stress 0.2689968 
Run 79 stress 0.2553429 
... Procrustes: rmse 0.003578207  max resid 0.02095359 
Run 80 stress 0.2624784 
*** Best solution was not repeated -- monoMDS stopping criteria:
    12: no. of iterations >= maxit
    68: stress ratio > sratmax

stressplot(df1s.mds2, main="Shepard plot")

df1s.mds2

Call:
metaMDS(comm = df1s.bc, k = 2, try = 40, trymax = 80, autotransform = FALSE,      maxit = 200) 

global Multidimensional Scaling using monoMDS

Data:     df1s.bc 
Distance: bray 

Dimensions: 2 
Stress:     0.2552758 
Stress type 1, weak ties
Best solution was not repeated after 80 tries
The best solution was from try 23 (random start)
Scaling: centring, PC rotation, halfchange scaling 
Species: scores missing

Stress values suggest 3d plot is a more accurate representation

Need to plot dim 1 vs dim 2 and dim 1 vs dim 3 with tsf groups identified by symbols

Generate graphs

This is a variant of code used in worked example for Ch 16. It does use the graph settings defined in appearance.R, which is called earlier.

a<-as.data.frame(scores(df1s.mds3))
a<-cbind(df[c(1:2)],a)   #Add site names & symbols from original data file
p1<-ggplot(data=a, aes(x=NMDS1, y=NMDS2, color=tsf, ) )+
  geom_point()+
  labs(y="MDS2", x="MDS1")+
  scale_shape_manual(values=sym3,
                     name="TSF",
                     guide =
                         guide_legend(label.theme = element_text(size=6), 
                                      title=NULL)
                     )
p2<-ggplot(data=a, aes(x=NMDS1, y=NMDS3, color = tsf, ) )+
  geom_point()+
  labs(y="MDS3", x="MDS1")+
  scale_shape_manual(values=sym3,
                     name="TSF",
                     guide =
                         guide_legend(label.theme = element_text(size=6),
                                    title=NULL)
                     )

p1 + p2 + plot_layout(guides='collect') & theme_qk() & scale_color_uchicago()

How would you fit a linear model to assess fire effects?

Do permanova including tsf

df1s.ado <- adonis2(df1s.bc~tsf,data=df,permutations=999)
print(df1s.ado)

Permutation test for adonis under reduced model
Terms added sequentially (first to last)
Permutation: free
Number of permutations: 999

adonis2(formula = df1s.bc ~ tsf, data = df, permutations = 999)
         Df SumOfSqs      R2      F Pr(>F)    
tsf       2   2.6433 0.11039 4.7151  0.001 ***
Residual 76  21.3028 0.88961                  
Total    78  23.9460 1.00000                  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Permanova suggests strong signal from tsf

SIMPER analysis for tsf group differences

df.sim <- simper(df1s[,-(1:2)], df$tsf, permutations=1000)
summary(df.sim)

Contrast: long_short 

     average      sd   ratio     ava     avb cumsum      p    
lgui 0.13247 0.14056 0.94250 0.15918 0.26922  0.172 0.3407    
eul  0.09332 0.12167 0.76700 0.02379 0.16793  0.293 0.0519 .  
dcor 0.08201 0.11054 0.74190 0.06582 0.10672  0.399 0.1279    
amac 0.06242 0.05874 1.06250 0.07712 0.07023  0.480 0.9610    
tnig 0.05948 0.07551 0.78760 0.03804 0.07852  0.557 0.8661    
ldel 0.05168 0.04157 1.24340 0.08687 0.00000  0.624 0.0010 ***
pent 0.05090 0.04718 1.07890 0.10975 0.03392  0.690 0.7373    
htal 0.03859 0.05635 0.68490 0.04599 0.02950  0.740 0.4006    
lwhi 0.03208 0.03938 0.81470 0.04549 0.01261  0.782 0.0230 *  
pspe 0.02827 0.03740 0.75580 0.04834 0.00680  0.819 0.0649 .  
amur 0.02647 0.05251 0.50400 0.04512 0.00000  0.853 0.8941    
rdie 0.02561 0.09803 0.26120 0.00512 0.04275  0.886 0.5125    
ecun 0.02190 0.04456 0.49150 0.03603 0.00000  0.915 0.0460 *  
esax 0.01860 0.04544 0.40940 0.01674 0.01687  0.939 0.1059    
ptex 0.01585 0.04231 0.37450 0.01170 0.01687  0.959 0.2298    
aram 0.01495 0.03676 0.40680 0.01517 0.01261  0.979 0.5764    
vros 0.01279 0.03799 0.33650 0.02216 0.00000  0.995 0.8412    
ppor 0.00363 0.01891 0.19220 0.00620 0.00000  1.000 0.7113    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Contrast: long_medium 

     average      sd   ratio     ava     avb cumsum      p    
lgui 0.10221 0.10333 0.98920 0.15918 0.21705  0.133 0.9820    
amac 0.10139 0.11774 0.86110 0.07712 0.16122  0.264 0.2138    
tnig 0.07029 0.10274 0.68410 0.03804 0.10788  0.356 0.6583    
pent 0.06601 0.06862 0.96190 0.10975 0.08112  0.442 0.0140 *  
amur 0.05952 0.11624 0.51210 0.04512 0.07283  0.519 0.0829 .  
dcor 0.05521 0.09128 0.60480 0.06582 0.03801  0.591 0.7842    
ldel 0.04942 0.03935 1.25580 0.08687 0.00185  0.655 0.0010 ***
eul  0.04752 0.07141 0.66540 0.02379 0.07484  0.716 0.9970    
htal 0.03907 0.06803 0.57430 0.04599 0.02956  0.767 0.3267    
vros 0.03267 0.09018 0.36230 0.02216 0.03599  0.809 0.2198    
lwhi 0.02992 0.03771 0.79360 0.04549 0.00840  0.848 0.0010 ***
pspe 0.02952 0.03909 0.75530 0.04834 0.00806  0.887 0.0020 ** 
ecun 0.02116 0.04260 0.49670 0.03603 0.00000  0.914 0.0010 ***
rdie 0.01844 0.06391 0.28850 0.00512 0.02970  0.938 0.8092    
aram 0.01711 0.04409 0.38800 0.01517 0.01739  0.960 0.5025    
esax 0.01138 0.02737 0.41600 0.01674 0.00420  0.975 0.6743    
ptex 0.00983 0.02516 0.39070 0.01170 0.00658  0.988 0.7193    
ppor 0.00929 0.03741 0.24840 0.00620 0.01107  1.000 0.3666    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Contrast: short_medium 

     average      sd   ratio     ava     avb cumsum      p   
lgui 0.15998 0.16074 0.99530 0.26922 0.21705  0.205 0.0060 **
amac 0.10618 0.14346 0.74010 0.07023 0.16122  0.340 0.1628   
eul  0.10520 0.12391 0.84900 0.16793 0.07484  0.475 0.0030 **
tnig 0.08903 0.11469 0.77630 0.07852 0.10788  0.589 0.0819 . 
dcor 0.07378 0.14450 0.51060 0.10672 0.03801  0.683 0.2378   
pent 0.04966 0.08217 0.60440 0.03392 0.08112  0.747 0.7992   
amur 0.04038 0.12310 0.32800 0.00000 0.07283  0.798 0.6883   
rdie 0.03751 0.11223 0.33420 0.04275 0.02970  0.847 0.1808   
htal 0.03122 0.07878 0.39630 0.02950 0.02956  0.886 0.7253   
vros 0.02238 0.09233 0.24240 0.00000 0.03599  0.915 0.6114   
aram 0.01559 0.04775 0.32650 0.01261 0.01739  0.935 0.5485   
ptex 0.01278 0.04462 0.28650 0.01687 0.00658  0.951 0.4156   
esax 0.01161 0.04217 0.27530 0.01687 0.00420  0.966 0.5754   
lwhi 0.01095 0.03778 0.28970 0.01261 0.00840  0.980 0.9930   
pspe 0.00840 0.02396 0.35050 0.00680 0.00806  0.991 0.9980   
ppor 0.00608 0.03495 0.17380 0.00000 0.01107  0.999 0.5495   
ldel 0.00101 0.00583 0.17380 0.00000 0.00185  1.000 1.0000   
ecun 0.00000 0.00000     NaN 0.00000 0.00000  1.000 0.9990   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Permutation: free
Number of permutations: 1000

Suppose you think of a reptile assemblage as simply the species that are present, regardless of how common they are. What would you conclude about the infuence of time since fire now?

Repeat the preceding analysis after creating a new data file that is presence-absence only.

df1.jac <- vegdist(df1,binary=TRUE,method='jaccard')
#Run 3d first
df1s.mds3 <- metaMDS(df1.jac,k=3,autotransform=FALSE,try=40,trymax=80,maxit=200)

Run 0 stress 0.1404128 
Run 1 stress 0.1480196 
Run 2 stress 0.1405858 
... Procrustes: rmse 0.01565156  max resid 0.08607928 
Run 3 stress 0.1485349 
Run 4 stress 0.1417498 
Run 5 stress 0.1404104 
... New best solution
... Procrustes: rmse 0.0023304  max resid 0.01808758 
Run 6 stress 0.1404131 
... Procrustes: rmse 0.002334262  max resid 0.01834106 
Run 7 stress 0.1415289 
Run 8 stress 0.140411 
... Procrustes: rmse 0.0007835425  max resid 0.005200482 
... Similar to previous best
Run 9 stress 0.1486406 
Run 10 stress 0.1404103 
... New best solution
... Procrustes: rmse 0.001207819  max resid 0.009469643 
... Similar to previous best
Run 11 stress 0.1404117 
... Procrustes: rmse 0.000916645  max resid 0.006636852 
... Similar to previous best
Run 12 stress 0.1405581 
... Procrustes: rmse 0.01501831  max resid 0.08924602 
Run 13 stress 0.143846 
Run 14 stress 0.1495992 
Run 15 stress 0.140558 
... Procrustes: rmse 0.01458579  max resid 0.08787773 
Run 16 stress 0.1404148 
... Procrustes: rmse 0.0015373  max resid 0.01187026 
Run 17 stress 0.1448137 
Run 18 stress 0.1443774 
Run 19 stress 0.140411 
... Procrustes: rmse 0.001777414  max resid 0.01404256 
Run 20 stress 0.1500017 
Run 21 stress 0.1485371 
Run 22 stress 0.1454737 
Run 23 stress 0.140411 
... Procrustes: rmse 0.00189255  max resid 0.01470197 
Run 24 stress 0.1496201 
Run 25 stress 0.1405615 
... Procrustes: rmse 0.01445625  max resid 0.08690181 
Run 26 stress 0.1477962 
Run 27 stress 0.1478653 
Run 28 stress 0.1488089 
Run 29 stress 0.1444777 
Run 30 stress 0.1405647 
... Procrustes: rmse 0.01491986  max resid 0.08994996 
Run 31 stress 0.1405594 
... Procrustes: rmse 0.01486184  max resid 0.08799154 
Run 32 stress 0.1405612 
... Procrustes: rmse 0.01487556  max resid 0.0876917 
Run 33 stress 0.140559 
... Procrustes: rmse 0.01485206  max resid 0.0880716 
Run 34 stress 0.1425376 
Run 35 stress 0.140417 
... Procrustes: rmse 0.001866647  max resid 0.01406319 
Run 36 stress 0.1418405 
Run 37 stress 0.1422689 
Run 38 stress 0.1486242 
Run 39 stress 0.1413517 
Run 40 stress 0.1468941 
*** Best solution repeated 2 times

stressplot(df1s.mds3, main="Shepard plot")

df1s.mds3

Call:
metaMDS(comm = df1.jac, k = 3, try = 40, trymax = 80, autotransform = FALSE,      maxit = 200) 

global Multidimensional Scaling using monoMDS

Data:     df1.jac 
Distance: binary jaccard 

Dimensions: 3 
Stress:     0.1404103 
Stress type 1, weak ties
Best solution was repeated 2 times in 40 tries
The best solution was from try 10 (random start)
Scaling: centring, PC rotation, halfchange scaling 
Species: scores missing

# Now look at 2d
df1s.mds2 <- metaMDS(df1.jac,k=2,autotransform=FALSE,try=40,trymax=80,maxit=200)

Run 0 stress 0.2009509 
Run 1 stress 0.2079375 
Run 2 stress 0.2224574 
Run 3 stress 0.2162073 
Run 4 stress 0.225869 
Run 5 stress 0.2284785 
Run 6 stress 0.2305984 
Run 7 stress 0.2199417 
Run 8 stress 0.202096 
Run 9 stress 0.2011626 
... Procrustes: rmse 0.006551419  max resid 0.03916879 
Run 10 stress 0.2104293 
Run 11 stress 0.2015109 
Run 12 stress 0.2009556 
... Procrustes: rmse 0.001703517  max resid 0.01335007 
Run 13 stress 0.2135392 
Run 14 stress 0.2246277 
Run 15 stress 0.2240908 
Run 16 stress 0.2250626 
Run 17 stress 0.200954 
... Procrustes: rmse 0.002998174  max resid 0.01914574 
Run 18 stress 0.2095741 
Run 19 stress 0.2186864 
Run 20 stress 0.203307 
Run 21 stress 0.2081138 
Run 22 stress 0.2212435 
Run 23 stress 0.2132584 
Run 24 stress 0.2055683 
Run 25 stress 0.2019164 
Run 26 stress 0.2108398 
Run 27 stress 0.2151842 
Run 28 stress 0.2009453 
... New best solution
... Procrustes: rmse 0.002748903  max resid 0.0197823 
Run 29 stress 0.206872 
Run 30 stress 0.2224539 
Run 31 stress 0.2031111 
Run 32 stress 0.2064081 
Run 33 stress 0.2149136 
Run 34 stress 0.2009451 
... New best solution
... Procrustes: rmse 0.002076934  max resid 0.0168433 
Run 35 stress 0.2111817 
Run 36 stress 0.2164562 
Run 37 stress 0.2027405 
Run 38 stress 0.2161414 
Run 39 stress 0.2242858 
Run 40 stress 0.2071548 
Run 41 stress 0.203393 
Run 42 stress 0.2266633 
Run 43 stress 0.2121352 
Run 44 stress 0.2175283 
Run 45 stress 0.2104325 
Run 46 stress 0.2029998 
Run 47 stress 0.2159185 
Run 48 stress 0.2042398 
Run 49 stress 0.2134569 
Run 50 stress 0.2145411 
Run 51 stress 0.2081612 
Run 52 stress 0.223412 
Run 53 stress 0.2327611 
Run 54 stress 0.2199861 
Run 55 stress 0.2115448 
Run 56 stress 0.209779 
Run 57 stress 0.2019152 
Run 58 stress 0.2236233 
Run 59 stress 0.2015106 
Run 60 stress 0.2036071 
Run 61 stress 0.2118656 
Run 62 stress 0.2120415 
Run 63 stress 0.2222671 
Run 64 stress 0.2176459 
Run 65 stress 0.2133805 
Run 66 stress 0.2009556 
... Procrustes: rmse 0.003462696  max resid 0.0201114 
Run 67 stress 0.201907 
Run 68 stress 0.2009456 
... Procrustes: rmse 0.0003121967  max resid 0.001856644 
... Similar to previous best
*** Best solution repeated 1 times

stressplot(df1s.mds2, main="Shepard plot")

df1s.mds2

Call:
metaMDS(comm = df1.jac, k = 2, try = 40, trymax = 80, autotransform = FALSE,      maxit = 200) 

global Multidimensional Scaling using monoMDS

Data:     df1.jac 
Distance: binary jaccard 

Dimensions: 2 
Stress:     0.2009451 
Stress type 1, weak ties
Best solution was repeated 1 time in 68 tries
The best solution was from try 34 (random start)
Scaling: centring, PC rotation, halfchange scaling 
Species: scores missing

Marginal decision here - stress just over 0.2 for k=2, so might be OK. Stress good for k=3, but 2-d MDS usually better for reporting or showing to audience

a<-as.data.frame(scores(df1s.mds3))
a<-cbind(df[c(1:2)],a)   #Add site names & symbols from original data file
p1<-ggplot(data=a, aes(x=NMDS1, y=NMDS2, color=tsf, ) )+
  geom_point()+
  labs(y="MDS2", x="MDS1")+
  scale_shape_manual(values=sym3,
                     name="TSF",
                     guide =
                         guide_legend(label.theme = element_text(size=6), 
                                      title=NULL)
                     )
p2<-ggplot(data=a, aes(x=NMDS1, y=NMDS3, color = tsf, ) )+
  geom_point()+
  labs(y="MDS3", x="MDS1")+
  scale_shape_manual(values=sym3,
                     name="TSF",
                     guide =
                         guide_legend(label.theme = element_text(size=6),
                                    title=NULL)
                     )

p1 + p2 + plot_layout(guides='collect') & theme_qk() & scale_color_uchicago()

Unburned sites separate very clearly on MDS1. Pattern clearer than with abundances included

Do permanova including tsf

df1.ado <- adonis2(df1.jac~tsf,data=df,permutations=999)
print(df1.ado)

Permutation test for adonis under reduced model
Terms added sequentially (first to last)
Permutation: free
Number of permutations: 999

adonis2(formula = df1.jac ~ tsf, data = df, permutations = 999)
         Df SumOfSqs      R2      F Pr(>F)    
tsf       2   3.7853 0.20678 9.9057  0.001 ***
Residual 76  14.5210 0.79322                  
Total    78  18.3062 1.00000                  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Clear fire signal

Extension activity

Dixon and her colleagues also recorded a range of habitat variables, which we’ve collected for you in dixonenv.csv

dixonenv <- read_csv("data/dixonenv.csv")

Rows: 79 Columns: 14── Column specification ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr  (4): site, tsf, veg, aspect
dbl (10): lcwd, shrcov, grcov, litcov, litdep, rocks, elev, warm, cold, twi
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

head(dixonenv, 10)

They recorded Coarse Woody Debris, which was long-transformed, litter cover (litcov), % cover of groundcover (grcov), and % cover of shrubs (shrcov), and rock cover. In the original paper, Dixon et al. were comfortable that these predictors weren’t correlated.

How is reptile assemblage related to time since fire and these five habitat variables

Check continuous predictors

boxplot(dixonenv$lcwd)

boxplot(dixonenv$litcov)

boxplot(dixonenv$grcov)

boxplot(dixonenv$shrcov)

boxplot(dixonenv$rocks)

#transform rocks to log+1
dixonenv$lrocks <- log10(dixonenv$rocks+1)
# check homog of dispersions on double standardized abundances
dixonbiota1s.disp <- betadisper(dixonbiota1s.bc,dixonbiota$tsf)
anova(dixonbiota1s.disp)

Analysis of Variance Table

Response: Distances
          Df  Sum Sq  Mean Sq F value    Pr(>F)    
Groups     2 0.48563 0.242815  32.165 7.573e-11 ***
Residuals 76 0.57374 0.007549                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

do permanova including tsf plus continuous predictors

dixonbiota1s.ado <- adonis2(dixonbiota1s.bc~tsf+lcwd+litcov+grcov+shrcov+lrocks,data=dixonenv,permutations=999)
print(dixonbiota1s.ado)

Permutation test for adonis under reduced model
Terms added sequentially (first to last)
Permutation: free
Number of permutations: 999

adonis2(formula = dixonbiota1s.bc ~ tsf + lcwd + litcov + grcov + shrcov + lrocks, data = dixonenv, permutations = 999)
         Df SumOfSqs      R2      F Pr(>F)    
tsf       2   2.6433 0.11039 4.8306  0.001 ***
lcwd      1   0.4865 0.02032 1.7783  0.052 .  
litcov    1   0.4037 0.01686 1.4757  0.135    
grcov     1   0.3281 0.01370 1.1991  0.296    
shrcov    1   0.3291 0.01374 1.2028  0.277    
lrocks    1   0.3300 0.01378 1.2062  0.269    
Residual 71  19.4253 0.81121                  
Total    78  23.9460 1.00000                  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#Run analysis for presence-absence as well
dixonbiota1.jac.ado <- adonis2(dixonbiota1.jac~tsf+lcwd+litcov+grcov+shrcov+lrocks,data=dixonenv,permutations=999)
print(dixonbiota1.jac.ado)

Permutation test for adonis under reduced model
Terms added sequentially (first to last)
Permutation: free
Number of permutations: 999

adonis2(formula = dixonbiota1.jac ~ tsf + lcwd + litcov + grcov + shrcov + lrocks, data = dixonenv, permutations = 999)
         Df SumOfSqs      R2       F Pr(>F)    
tsf       2   3.7853 0.20678 10.4082  0.001 ***
lcwd      1   0.3802 0.02077  2.0908  0.043 *  
litcov    1   0.3145 0.01718  1.7296  0.097 .  
grcov     1   0.4097 0.02238  2.2530  0.021 *  
shrcov    1   0.2748 0.01501  1.5115  0.124    
lrocks    1   0.2309 0.01262  1.2700  0.233    
Residual 71  12.9108 0.70527                   
Total    78  18.3062 1.00000                   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Main effect is tsf; only some question about CWD when we look at species composition and abundance, while for presence-absence data, there are also effects of CWD and groundcover

mvabund

dixonbiotamv <- mvabund(dixonbiota[,-(1:2)])
dixonbiotamv.mv <- manyglm(dixonbiotamv~dixonenv$tsf+dixonenv$lcwd+dixonenv$litcov+dixonenv$grcov+dixonenv$shrcov+dixonenv$lrocks,family="poisson")
plot(dixonbiotamv.mv)

# still hetero vars so use -ve binomial
dixonbiotamv1.mv <- manyglm(dixonbiotamv~dixonenv$tsf+dixonenv$lcwd+dixonenv$litcov+dixonenv$grcov+dixonenv$shrcov+dixonenv$lrocks,family="negative_binomial")
plot(dixonbiotamv1.mv)

anova(dixonbiotamv1.mv)

Time elapsed: 0 hr 1 min 39 sec
Analysis of Deviance Table

Model: dixonbiotamv ~ dixonenv$tsf + dixonenv$lcwd + dixonenv$litcov + dixonenv$grcov + dixonenv$shrcov + dixonenv$lrocks

Multivariate test:
                Res.Df Df.diff   Dev Pr(>Dev)    
(Intercept)         78                           
dixonenv$tsf        76       2 331.6    0.001 ***
dixonenv$lcwd       75       1  61.5    0.001 ***
dixonenv$litcov     74       1  25.9    0.265    
dixonenv$grcov      73       1  36.4    0.046 *  
dixonenv$shrcov     72       1  21.5    0.607    
dixonenv$lrocks     71       1  58.0    0.001 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Arguments:
 Test statistics calculated assuming uncorrelated response (for faster computation) 
 P-value calculated using 999 iterations via PIT-trap resampling.

B Simple MDS & Permanova

Let’s return to the Hutto and Barrett (2021) example from the Chapter 15 exercises. There, we used similarity-based analyses to explore the relationship between frog assemblages in ponds and the degree of urbanization surrounding. This seems a question that could just as easily be examined using distances or dissimilarities.

Return to the data and assess whether frog assemblages (using the standardized abundance scale or presence-absence) differ with urbanization.

Did your conclusions differ from those obtained using RDA?

df <- read_csv("data/huttoamph.csv")

Rows: 50 Columns: 14── Column specification ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr  (1): type
dbl (13): Site, ACRE, BAME, BFOW, GCAR, HCIN, HVER, LCAT, LCLA, LPAL, LSPH, PCRU, PFER
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

head(df,10)

#ponds 9 and 40 had no frogs. Remove for distance analysis
df <- df %>% 
  filter(Site!= 9 & Site !=40)
df

The explanation of the variables is in the previous chapter’s exercises.

df1 <- df[,-(1:2)]   #Stripped back file without labels, leaving only 
df1.bc <- vegdist(df1,'bray')

#Run 2d first
df1.mds2 <- metaMDS(df1.bc,k=2,autotransform=FALSE,try=40,trymax=80,maxit=200)

Run 0 stress 0.196708 
Run 1 stress 0.1877626 
... New best solution
... Procrustes: rmse 0.06205021  max resid 0.2289634 
Run 2 stress 0.1888251 
Run 3 stress 0.1920872 
Run 4 stress 0.1853642 
... New best solution
... Procrustes: rmse 0.07936891  max resid 0.4999309 
Run 5 stress 0.1855846 
... Procrustes: rmse 0.01953946  max resid 0.1056279 
Run 6 stress 0.1984394 
Run 7 stress 0.1862238 
Run 8 stress 0.200767 
Run 9 stress 0.1961235 
Run 10 stress 0.1853641 
... New best solution
... Procrustes: rmse 9.823368e-05  max resid 0.0004897303 
... Similar to previous best
Run 11 stress 0.204919 
Run 12 stress 0.2001183 
Run 13 stress 0.1875595 
Run 14 stress 0.2031338 
Run 15 stress 0.1875596 
Run 16 stress 0.2001202 
Run 17 stress 0.1877626 
Run 18 stress 0.1968146 
Run 19 stress 0.2191092 
Run 20 stress 0.1862237 
Run 21 stress 0.1862237 
Run 22 stress 0.2016092 
Run 23 stress 0.1919379 
Run 24 stress 0.1921465 
Run 25 stress 0.1875596 
Run 26 stress 0.1862241 
Run 27 stress 0.1888273 
Run 28 stress 0.1958875 
Run 29 stress 0.2103827 
Run 30 stress 0.196403 
Run 31 stress 0.1920872 
Run 32 stress 0.1853644 
... Procrustes: rmse 0.0001814753  max resid 0.0007769668 
... Similar to previous best
Run 33 stress 0.1929083 
Run 34 stress 0.1948815 
Run 35 stress 0.2023845 
Run 36 stress 0.1961143 
Run 37 stress 0.1978435 
Run 38 stress 0.214542 
Run 39 stress 0.1970125 
Run 40 stress 0.2055679 
*** Best solution repeated 2 times

stressplot(df1.mds2, main="Shepard plot")

df1.mds2

Call:
metaMDS(comm = df1.bc, k = 2, try = 40, trymax = 80, autotransform = FALSE,      maxit = 200) 

global Multidimensional Scaling using monoMDS

Data:     df1.bc 
Distance: bray 

Dimensions: 2 
Stress:     0.1853641 
Stress type 1, weak ties
Best solution was repeated 2 times in 40 tries
The best solution was from try 10 (random start)
Scaling: centring, PC rotation, halfchange scaling 
Species: scores missing

# Now look at 3d
df.mds3 <- metaMDS(df1.bc,k=3,autotransform=FALSE,try=40,trymax=80,maxit=200)

Run 0 stress 0.129187 
Run 1 stress 0.1388768 
Run 2 stress 0.1382418 
Run 3 stress 0.1383885 
Run 4 stress 0.129186 
... New best solution
... Procrustes: rmse 0.001355541  max resid 0.007232405 
... Similar to previous best
Run 5 stress 0.129186 
... Procrustes: rmse 0.0005937864  max resid 0.001870774 
... Similar to previous best
Run 6 stress 0.1291858 
... New best solution
... Procrustes: rmse 0.0008376747  max resid 0.004337427 
... Similar to previous best
Run 7 stress 0.129186 
... Procrustes: rmse 0.0007498725  max resid 0.00420023 
... Similar to previous best
Run 8 stress 0.1291864 
... Procrustes: rmse 0.0002763668  max resid 0.0009518002 
... Similar to previous best
Run 9 stress 0.1383896 
Run 10 stress 0.1384317 
Run 11 stress 0.1289614 
... New best solution
... Procrustes: rmse 0.0109308  max resid 0.05605748 
Run 12 stress 0.128962 
... Procrustes: rmse 0.0009039712  max resid 0.005251417 
... Similar to previous best
Run 13 stress 0.1383929 
Run 14 stress 0.1402552 
Run 15 stress 0.128962 
... Procrustes: rmse 0.0009246251  max resid 0.005329771 
... Similar to previous best
Run 16 stress 0.1291859 
... Procrustes: rmse 0.01058901  max resid 0.0547769 
Run 17 stress 0.1289615 
... Procrustes: rmse 0.000419003  max resid 0.002034791 
... Similar to previous best
Run 18 stress 0.1380889 
Run 19 stress 0.129186 
... Procrustes: rmse 0.01092673  max resid 0.05630155 
Run 20 stress 0.1383886 
Run 21 stress 0.1289617 
... Procrustes: rmse 0.0002465603  max resid 0.0007331739 
... Similar to previous best
Run 22 stress 0.1291604 
... Procrustes: rmse 0.009890855  max resid 0.05184425 
Run 23 stress 0.1291601 
... Procrustes: rmse 0.009490803  max resid 0.05078321 
Run 24 stress 0.1291855 
... Procrustes: rmse 0.01066625  max resid 0.05496251 
Run 25 stress 0.1391225 
Run 26 stress 0.129186 
... Procrustes: rmse 0.01099671  max resid 0.05629543 
Run 27 stress 0.1291866 
... Procrustes: rmse 0.01114331  max resid 0.05711187 
Run 28 stress 0.1384282 
Run 29 stress 0.1291861 
... Procrustes: rmse 0.01101164  max resid 0.0563452 
Run 30 stress 0.1383979 
Run 31 stress 0.1289613 
... New best solution
... Procrustes: rmse 0.0004223583  max resid 0.001654659 
... Similar to previous best
Run 32 stress 0.1380874 
Run 33 stress 0.1291858 
... Procrustes: rmse 0.01077276  max resid 0.05558138 
Run 34 stress 0.1291859 
... Procrustes: rmse 0.01081878  max resid 0.05582879 
Run 35 stress 0.1291858 
... Procrustes: rmse 0.01042548  max resid 0.05433355 
Run 36 stress 0.1385874 
Run 37 stress 0.1291863 
... Procrustes: rmse 0.01092361  max resid 0.05638233 
Run 38 stress 0.1308828 
Run 39 stress 0.1384275 
Run 40 stress 0.1388191 
*** Best solution repeated 1 times

stressplot(df.mds3, main="Shepard plot")

df.mds3

Call:
metaMDS(comm = df1.bc, k = 3, try = 40, trymax = 80, autotransform = FALSE,      maxit = 200) 

global Multidimensional Scaling using monoMDS

Data:     df1.bc 
Distance: bray 

Dimensions: 3 
Stress:     0.1289613 
Stress type 1, weak ties
Best solution was repeated 1 time in 40 tries
The best solution was from try 31 (random start)
Scaling: centring, PC rotation, halfchange scaling 
Species: scores missing

2-d is acceptable

a<-as.data.frame(scores(df1.mds2))
a<-cbind(df[c(1:2)],a)   #Add site names & symbols from original data file
ggplot(data=a, aes(x=NMDS1, y=NMDS2, color=type, ) )+
  geom_point()+
  labs(y="MDS2", x="MDS1")+
  scale_shape_manual(values=sym3,
                     name="Urbanization",
                     guide =
                         guide_legend(label.theme = element_text(size=6), 
                                      title=NULL)
                     )+
  theme_qk() + scale_color_uchicago()

df1.ado <- adonis2(df1.bc~type,data=df,permutations=999)
print(df1.ado)

Permutation test for adonis under reduced model
Terms added sequentially (first to last)
Permutation: free
Number of permutations: 999

adonis2(formula = df1.bc ~ type, data = df, permutations = 999)
         Df SumOfSqs      R2      F Pr(>F)  
type      2    0.912 0.07281 1.7669  0.045 *
Residual 45   11.614 0.92719                
Total    47   12.526 1.00000                
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Permanova detects a difference; MDS plot suggests low urbanization may be different from the other two, though worth looking at dispersion as well, as spread of this type is less than the other two on plot

Repeat analysis with presence-absence

df1.jac <- vegdist(df1,binary=TRUE,method=‘jaccard’)

#Run 2d first
df1.jac <- vegdist(df1,binary=TRUE,method='jaccard')
df1.mds2 <- metaMDS(df1.jac,k=2,autotransform=FALSE,try=40,trymax=80,maxit=200)

Run 0 stress 0.1551612 
Run 1 stress 0.1630921 
Run 2 stress 0.1622263 
Run 3 stress 0.1724054 
Run 4 stress 0.1585799 
Run 5 stress 0.1761433 
Run 6 stress 0.1608258 
Run 7 stress 0.1734443 
Run 8 stress 0.1621126 
Run 9 stress 0.1536788 
... New best solution
... Procrustes: rmse 0.01901219  max resid 0.06506965 
Run 10 stress 0.1760283 
Run 11 stress 0.1609339 
Run 12 stress 0.1565712 
Run 13 stress 0.1939366 
Run 14 stress 0.1843833 
Run 15 stress 0.162526 
Run 16 stress 0.1622263 
Run 17 stress 0.1570453 
Run 18 stress 0.1708675 
Run 19 stress 0.1537515 
... Procrustes: rmse 0.003076067  max resid 0.01474863 
Run 20 stress 0.1539129 
... Procrustes: rmse 0.00947995  max resid 0.04425278 
Run 21 stress 0.1664178 
Run 22 stress 0.1795028 
Run 23 stress 0.1585801 
Run 24 stress 0.1555104 
Run 25 stress 0.1708671 
Run 26 stress 0.167874 
Run 27 stress 0.1673092 
Run 28 stress 0.1719174 
Run 29 stress 0.1591019 
Run 30 stress 0.1711482 
Run 31 stress 0.1795847 
Run 32 stress 0.1608763 
Run 33 stress 0.1721431 
Run 34 stress 0.1585798 
Run 35 stress 0.1850134 
Run 36 stress 0.1607836 
Run 37 stress 0.1586731 
Run 38 stress 0.1738412 
Run 39 stress 0.2004114 
Run 40 stress 0.1672269 
Run 41 stress 0.1585799 
Run 42 stress 0.1591019 
Run 43 stress 0.162564 
Run 44 stress 0.1620556 
Run 45 stress 0.1679859 
Run 46 stress 0.171033 
Run 47 stress 0.1617404 
Run 48 stress 0.1538246 
... Procrustes: rmse 0.01103138  max resid 0.05221548 
Run 49 stress 0.172142 
Run 50 stress 0.1960733 
Run 51 stress 0.1731862 
Run 52 stress 0.1627813 
Run 53 stress 0.1536788 
... Procrustes: rmse 0.0001144048  max resid 0.000322326 
... Similar to previous best
*** Best solution repeated 1 times

stressplot(df1.mds2, main="Shepard plot")

df1.mds2

Call:
metaMDS(comm = df1.jac, k = 2, try = 40, trymax = 80, autotransform = FALSE,      maxit = 200) 

global Multidimensional Scaling using monoMDS

Data:     df1.jac 
Distance: binary jaccard 

Dimensions: 2 
Stress:     0.1536788 
Stress type 1, weak ties
Best solution was repeated 1 time in 53 tries
The best solution was from try 9 (random start)
Scaling: centring, PC rotation, halfchange scaling 
Species: scores missing

# Now look at 3d
df.mds3 <- metaMDS(df1.jac,k=3,autotransform=FALSE,try=40,trymax=80,maxit=200)

Run 0 stress 0.09466813 
Run 1 stress 0.09466802 
... New best solution
... Procrustes: rmse 0.0004764601  max resid 0.002102695 
... Similar to previous best
Run 2 stress 0.09466791 
... New best solution
... Procrustes: rmse 0.0001518385  max resid 0.0004555261 
... Similar to previous best
Run 3 stress 0.09466803 
... Procrustes: rmse 0.0001108676  max resid 0.0004745956 
... Similar to previous best
Run 4 stress 0.1005052 
Run 5 stress 0.09466787 
... New best solution
... Procrustes: rmse 5.340853e-05  max resid 0.0001992126 
... Similar to previous best
Run 6 stress 0.1005044 
Run 7 stress 0.1005056 
Run 8 stress 0.09466803 
... Procrustes: rmse 0.0001288817  max resid 0.0005023949 
... Similar to previous best
Run 9 stress 0.1005047 
Run 10 stress 0.09466793 
... Procrustes: rmse 5.310415e-05  max resid 0.0002709636 
... Similar to previous best
Run 11 stress 0.09466796 
... Procrustes: rmse 0.0001256892  max resid 0.0003435652 
... Similar to previous best
Run 12 stress 0.094668 
... Procrustes: rmse 9.745181e-05  max resid 0.0004669281 
... Similar to previous best
Run 13 stress 0.09693819 
Run 14 stress 0.09694632 
Run 15 stress 0.1005051 
Run 16 stress 0.09466833 
... Procrustes: rmse 0.0002133226  max resid 0.001028399 
... Similar to previous best
Run 17 stress 0.09466821 
... Procrustes: rmse 0.0002054658  max resid 0.0009036094 
... Similar to previous best
Run 18 stress 0.09466811 
... Procrustes: rmse 0.0001423043  max resid 0.0004843394 
... Similar to previous best
Run 19 stress 0.09466805 
... Procrustes: rmse 0.0003941737  max resid 0.002017784 
... Similar to previous best
Run 20 stress 0.109714 
Run 21 stress 0.1005051 
Run 22 stress 0.09466816 
... Procrustes: rmse 0.000195071  max resid 0.0006552548 
... Similar to previous best
Run 23 stress 0.09466783 
... New best solution
... Procrustes: rmse 6.655458e-05  max resid 0.0002929959 
... Similar to previous best
Run 24 stress 0.09466805 
... Procrustes: rmse 0.0001689595  max resid 0.0006247545 
... Similar to previous best
Run 25 stress 0.1005048 
Run 26 stress 0.09466801 
... Procrustes: rmse 0.0001773767  max resid 0.0009243069 
... Similar to previous best
Run 27 stress 0.09694616 
Run 28 stress 0.1004898 
Run 29 stress 0.09466783 
... New best solution
... Procrustes: rmse 0.0001519669  max resid 0.0004171994 
... Similar to previous best
Run 30 stress 0.1005053 
Run 31 stress 0.1097146 
Run 32 stress 0.1004899 
Run 33 stress 0.09466786 
... Procrustes: rmse 0.0001828158  max resid 0.000484661 
... Similar to previous best
Run 34 stress 0.1090648 
Run 35 stress 0.09466812 
... Procrustes: rmse 0.0002666448  max resid 0.0009620816 
... Similar to previous best
Run 36 stress 0.09694549 
Run 37 stress 0.09694624 
Run 38 stress 0.094668 
... Procrustes: rmse 0.0002647032  max resid 0.001383475 
... Similar to previous best
Run 39 stress 0.09466798 
... Procrustes: rmse 0.0002442504  max resid 0.0007524015 
... Similar to previous best
Run 40 stress 0.09466802 
... Procrustes: rmse 0.000254564  max resid 0.0008115543 
... Similar to previous best
*** Best solution repeated 6 times

stressplot(df.mds3, main="Shepard plot")

df.mds3

Call:
metaMDS(comm = df1.jac, k = 3, try = 40, trymax = 80, autotransform = FALSE,      maxit = 200) 

global Multidimensional Scaling using monoMDS

Data:     df1.jac 
Distance: binary jaccard 

Dimensions: 3 
Stress:     0.09466783 
Stress type 1, weak ties
Best solution was repeated 6 times in 40 tries
The best solution was from try 29 (random start)
Scaling: centring, PC rotation, halfchange scaling 
Species: scores missing

Again, stress OK with 2D, good with 3D. Look at 2D first

a<-as.data.frame(scores(df1.mds2))
a<-cbind(df[c(1:2)],a)   #Add site names & symbols from original data file
ggplot(data=a, aes(x=NMDS1, y=NMDS2, color=type, ) )+
  geom_point()+
  labs(y="MDS2", x="MDS1")+
  scale_shape_manual(values=sym3,
                     name="Urbanization",
                     guide =
                         guide_legend(label.theme = element_text(size=6), 
                                      title=NULL)
                     )+
  theme_qk() + scale_color_uchicago()

df1.ado <- adonis2(df1.jac~type,data=df,permutations=999)
print(df1.ado)

Permutation test for adonis under reduced model
Terms added sequentially (first to last)
Permutation: free
Number of permutations: 999

adonis2(formula = df1.jac ~ type, data = df, permutations = 999)
         Df SumOfSqs      R2      F Pr(>F)
type      2   0.4525 0.04389 1.0328  0.403
Residual 45   9.8587 0.95611              
Total    47  10.3112 1.00000              

No separation seen on Presence-Absence

C

Griffen et al. (2019) examined the oral microbiome of humans, with a focus on the effects of HIV and AntiRetroviral therapy (ART) on this microbiome. They described the microbiome of 341 patients, who fell into one of two categories: HIV^-, and HIV+ with ART. They also matched their samples as far as possible for sex, along with other conditions that can influence oral microbiomes (Candida infection, current smoking). We’ll use their records for these other categories as well. There were more HIV+ than - subjects, but within these groups, approximately equal numbers of two sexes, two Candida infection status, and two smoking categories.

The bacteriome was recorded separately using 16S RNA sequencing, which overed just over 600 taxa.

The data are available from their paper, as two Excel files in the supplementary information. The metadata gives HIV status, and a raft of demographic information. We’ll just work with HIV, sex, Candida, and current smoking. A second sheet has the bacteriome. The code chunk below reads this file (from a Downloads folder - you’ll need to modify the file location). It also screens the file for any bacterial taxa that were not present in this particular study, which reduces the bacterial taxa to 599.

library(readxl)
#Read in metadata. Lots of information we don't need, so a couple of iterations to get down to a few columns we want - HIV status, candida, current smoking, gender
df <- read_excel("~/Downloads/41598_2019_55703_MOESM3_ESM.xlsx", 
     range = "A1:S342", na = "NA")
df<-select(df, HIV, candida, gender, smoking_current)
head(df)

#df1 is the bacterial data
df1 <- read_excel("~/Downloads/41598_2019_55703_MOESM2_ESM.xlsx")

New names:

head(df1,10)

df1 <- df1 %>%
  select(where( ~ is.numeric(.x) && sum(.x) != 0))  #Drops any bacterial taxon not recorded in this study (i.e. where it's all zeroes)
df1 <- df1[,-1] #remove col1, which is sample ID

Use a dissimilarity based approach to assess the combined effects of HIV-ART, sex, Candida infection, and current smoking on the bacteriome

Do these factors act independently on the bacteriome?

Which effects are largest (use some graphical methods)?

First steps: Think about standardization and which distance measure you’ll use. Bray-Curtis seems common for these kind of data, and and counts of different taxa vary widely.

df1s <- wisconsin(df1)   #Common to do some standardisation, and counts of different taxa vary widely
df1s.bc <- vegdist(df1s,'bray')   #With standardisation (wisconsin)

Run factorial model using permanova, based on B-C

df1.ado <- adonis2(df1s.bc~HIV*candida*gender*smoking_current,data=df,permutations=999)
print(df1.ado)

Permutation test for adonis under reduced model
Terms added sequentially (first to last)
Permutation: free
Number of permutations: 999

adonis2(formula = df1s.bc ~ HIV * candida * gender * smoking_current, data = df, permutations = 999)
                                    Df SumOfSqs      R2      F Pr(>F)    
HIV                                  1    1.369 0.01250 4.3776  0.001 ***
candida                              1    1.461 0.01333 4.6694  0.001 ***
gender                               1    0.492 0.00449 1.5737  0.009 ** 
smoking_current                      1    1.111 0.01014 3.5513  0.001 ***
HIV:candida                          1    0.318 0.00290 1.0151  0.401    
HIV:gender                           1    0.354 0.00324 1.1332  0.204    
candida:gender                       1    0.257 0.00234 0.8202  0.861    
HIV:smoking_current                  1    0.269 0.00245 0.8584  0.765    
candida:smoking_current              1    0.344 0.00314 1.0989  0.243    
gender:smoking_current               1    0.394 0.00360 1.2602  0.098 .  
HIV:candida:gender                   1    0.349 0.00319 1.1164  0.203    
HIV:candida:smoking_current          1    0.280 0.00255 0.8943  0.707    
HIV:gender:smoking_current           1    0.265 0.00242 0.8463  0.793    
candida:gender:smoking_current       1    0.254 0.00232 0.8110  0.841    
HIV:candida:gender:smoking_current   1    0.362 0.00330 1.1562  0.213    
Residual                           325  101.664 0.92809                  
Total                              340  109.541 1.00000                  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Model shows main effects, but no interactions important (or at least, detected)

Could you fit a simpler model to the data?

Run simpler model

df2.ado <- adonis2(df1s.bc~HIV+candida+gender+smoking_current,data=df,permutations=999)
print(df2.ado)

Permutation test for adonis under reduced model
Terms added sequentially (first to last)
Permutation: free
Number of permutations: 999

adonis2(formula = df1s.bc ~ HIV + candida + gender + smoking_current, data = df, permutations = 999)
                 Df SumOfSqs      R2      F Pr(>F)    
HIV               1    1.369 0.01250 4.3775  0.001 ***
candida           1    1.461 0.01333 4.6693  0.001 ***
gender            1    0.492 0.00449 1.5736  0.013 *  
smoking_current   1    1.111 0.01014 3.5512  0.001 ***
Residual        336  105.108 0.95953                  
Total           340  109.541 1.00000                  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

No change; not too surprising, as all four main effects were detected easily.

Now go to data visualization

#Run 3d first
df1s.mds3 <- metaMDS(df1s.bc,k=3,autotransform=FALSE,try=40,trymax=80,maxit=200)

Run 0 stress 0.1573956 
Run 1 stress 0.1579391 
Run 2 stress 0.1572758 
... New best solution
... Procrustes: rmse 0.01025549  max resid 0.1525195 
Run 3 stress 0.1579179 
Run 4 stress 0.157493 
... Procrustes: rmse 0.007283638  max resid 0.1009598 
Run 5 stress 0.1578733 
Run 6 stress 0.1593345 
Run 7 stress 0.1596534 
Run 8 stress 0.1578064 
Run 9 stress 0.1576408 
... Procrustes: rmse 0.005563033  max resid 0.08072257 
Run 10 stress 0.1574589 
... Procrustes: rmse 0.007053274  max resid 0.1007079 
Run 11 stress 0.1574113 
... Procrustes: rmse 0.008390411  max resid 0.1530435 
Run 12 stress 0.1584759 
Run 13 stress 0.1590917 
Run 14 stress 0.1579313 
Run 15 stress 0.1577686 
... Procrustes: rmse 0.00731071  max resid 0.08902012 
Run 16 stress 0.1579681 
Run 17 stress 0.157889 
Run 18 stress 0.1576122 
... Procrustes: rmse 0.01046473  max resid 0.1079021 
Run 19 stress 0.1575193 
... Procrustes: rmse 0.01344117  max resid 0.1517403 
Run 20 stress 0.1590183 
Run 21 stress 0.1572873 
... Procrustes: rmse 0.005938186  max resid 0.1084724 
Run 22 stress 0.1579056 
Run 23 stress 0.1578609 
Run 24 stress 0.1575202 
... Procrustes: rmse 0.01348549  max resid 0.1517058 
Run 25 stress 0.1580688 
Run 26 stress 0.1573537 
... Procrustes: rmse 0.005408552  max resid 0.09851134 
Run 27 stress 0.1584273 
Run 28 stress 0.1575528 
... Procrustes: rmse 0.007336861  max resid 0.1083012 
Run 29 stress 0.1591385 
Run 30 stress 0.1575529 
... Procrustes: rmse 0.007325558  max resid 0.1082856 
Run 31 stress 0.1584616 
Run 32 stress 0.157776 
Run 33 stress 0.1579148 
Run 34 stress 0.1590267 
Run 35 stress 0.1574802 
... Procrustes: rmse 0.01063831  max resid 0.1527981 
Run 36 stress 0.1574825 
... Procrustes: rmse 0.0125045  max resid 0.1520643 
Run 37 stress 0.158388 
Run 38 stress 0.1581063 
Run 39 stress 0.1586636 
Run 40 stress 0.159344 
Run 41 stress 0.1575585 
... Procrustes: rmse 0.007481761  max resid 0.1081395 
Run 42 stress 0.1580441 
Run 43 stress 0.1583809 
Run 44 stress 0.1586431 
Run 45 stress 0.1575179 
... Procrustes: rmse 0.009346751  max resid 0.1082067 
Run 46 stress 0.1585796 
Run 47 stress 0.1578686 
Run 48 stress 0.1595925 
Run 49 stress 0.1579636 
Run 50 stress 0.1572746 
... New best solution
... Procrustes: rmse 0.0005494158  max resid 0.00824335 
... Similar to previous best
*** Best solution repeated 1 times

stressplot(df1s.mds3, main="Shepard plot")

df1s.mds3

Call:
metaMDS(comm = df1s.bc, k = 3, try = 40, trymax = 80, autotransform = FALSE,      maxit = 200) 

global Multidimensional Scaling using monoMDS

Data:     df1s.bc 
Distance: bray 

Dimensions: 3 
Stress:     0.1572746 
Stress type 1, weak ties
Best solution was repeated 1 time in 50 tries
The best solution was from try 50 (random start)
Scaling: centring, PC rotation, halfchange scaling 
Species: scores missing

# Now look at 2d
df1s.mds2 <- metaMDS(df1s.bc,k=2,autotransform=FALSE,try=40,trymax=80,maxit=200)

Run 0 stress 0.2249552 
Run 1 stress 0.2254408 
... Procrustes: rmse 0.01613467  max resid 0.2308209 
Run 2 stress 0.2253412 
... Procrustes: rmse 0.01606138  max resid 0.2302704 
Run 3 stress 0.2250828 
... Procrustes: rmse 0.009539203  max resid 0.1512118 
Run 4 stress 0.2250839 
... Procrustes: rmse 0.009650024  max resid 0.1525654 
Run 5 stress 0.2250832 
... Procrustes: rmse 0.009612515  max resid 0.1518641 
Run 6 stress 0.2248224 
... New best solution
... Procrustes: rmse 0.00631171  max resid 0.08263297 
Run 7 stress 0.2253967 
Run 8 stress 0.225527 
Run 9 stress 0.2248128 
... New best solution
... Procrustes: rmse 0.0004834727  max resid 0.004762261 
... Similar to previous best
Run 10 stress 0.2255447 
Run 11 stress 0.2249136 
... Procrustes: rmse 0.00484646  max resid 0.08193293 
Run 12 stress 0.2500137 
Run 13 stress 0.2256121 
Run 14 stress 0.2254928 
Run 15 stress 0.2258491 
Run 16 stress 0.2248191 
... Procrustes: rmse 0.0008375058  max resid 0.008573354 
... Similar to previous best
Run 17 stress 0.225849 
Run 18 stress 0.225634 
Run 19 stress 0.2249543 
... Procrustes: rmse 0.006228313  max resid 0.08122835 
Run 20 stress 0.2250713 
... Procrustes: rmse 0.008311559  max resid 0.1485133 
Run 21 stress 0.2258052 
Run 22 stress 0.2248935 
... Procrustes: rmse 0.004390114  max resid 0.08026974 
Run 23 stress 0.2257226 
Run 24 stress 0.2253296 
Run 25 stress 0.2257815 
Run 26 stress 0.2254095 
Run 27 stress 0.2257822 
Run 28 stress 0.2248496 
... Procrustes: rmse 0.004509314  max resid 0.0818388 
Run 29 stress 0.2250698 
... Procrustes: rmse 0.008236929  max resid 0.1489688 
Run 30 stress 0.225152 
... Procrustes: rmse 0.01050374  max resid 0.1503719 
Run 31 stress 0.2426061 
Run 32 stress 0.2253598 
Run 33 stress 0.418577 
Run 34 stress 0.2254408 
Run 35 stress 0.2254011 
Run 36 stress 0.2257955 
Run 37 stress 0.2256214 
Run 38 stress 0.22512 
... Procrustes: rmse 0.009424238  max resid 0.1483398 
Run 39 stress 0.2255057 
Run 40 stress 0.225153 
... Procrustes: rmse 0.0106079  max resid 0.1526557 
*** Best solution repeated 2 times

stressplot(df1s.mds2, main="Shepard plot")

df1s.mds2

Call:
metaMDS(comm = df1s.bc, k = 2, try = 40, trymax = 80, autotransform = FALSE,      maxit = 200) 

global Multidimensional Scaling using monoMDS

Data:     df1s.bc 
Distance: bray 

Dimensions: 2 
Stress:     0.2248128 
Stress type 1, weak ties
Best solution was repeated 2 times in 40 tries
The best solution was from try 9 (random start)
Scaling: centring, PC rotation, halfchange scaling 
Species: scores missing

3D model fits better than 2D. Stress for 2D above 0.2

Start with plots where symbol colour indicates levels of one of the factors

a<-as.data.frame(scores(df1s.mds3))
a<-cbind(df[c(1:4)],a)   #Add site names & symbols from original data file
p1<-ggplot(data=a, aes(x=NMDS1, y=NMDS2, color=HIV) )+
  geom_point()+
  stat_ellipse()+
  labs(y="MDS2", x="MDS1")+
  scale_shape_manual(values=c(16,17),
                     guide =
                         guide_legend(label.theme = element_text(size=6), 
                                      title=NULL)
                     )
p2<-ggplot(data=a, aes(x=NMDS1, y=NMDS3, color = HIV) )+
  geom_point()+
  stat_ellipse()+
  labs(y="MDS3", x="MDS1")+
  scale_shape_manual(values=c(16,17),
                     guide =
                         guide_legend(label.theme = element_text(size=6),
                                    title=NULL)
                     )

p1 + p2 + plot_layout(guides='collect') + plot_annotation(title="HIV status") & theme_qk() & scale_color_uchicago()

p1<-ggplot(data=a, aes(x=NMDS1, y=NMDS2, color=candida) )+
  geom_point()+
  stat_ellipse()+
  labs(y="MDS2", x="MDS1")+
  scale_shape_manual(values=c(16,17),
                     guide =
                         guide_legend(label.theme = element_text(size=6), 
                                      title=NULL)
                     )
p2<-ggplot(data=a, aes(x=NMDS1, y=NMDS3, color = candida) )+
  geom_point()+
  stat_ellipse()+
  labs(y="MDS3", x="MDS1")+
  scale_shape_manual(values=c(16,17),
                     guide =
                         guide_legend(label.theme = element_text(size=6),
                                    title=NULL)
                     )

p1 + p2 + plot_layout(guides='collect') + plot_annotation(title="Candida") & theme_qk() & scale_color_uchicago()

p1<-ggplot(data=a, aes(x=NMDS1, y=NMDS2, color=gender) )+
  geom_point()+
  stat_ellipse()+
  labs(y="MDS2", x="MDS1")+
  scale_shape_manual(values=c(16,17),
                     guide =
                         guide_legend(label.theme = element_text(size=6), 
                                      title=NULL)
                     )
p2<-ggplot(data=a, aes(x=NMDS1, y=NMDS3, color = gender) )+
  geom_point()+
  stat_ellipse()+
  labs(y="MDS3", x="MDS1")+
  scale_shape_manual(values=c(16,17),
                     guide =
                         guide_legend(label.theme = element_text(size=6),
                                    title=NULL)
                     )

p1 + p2 + plot_layout(guides='collect') + plot_annotation(title="Sex") & theme_qk() & scale_color_uchicago()

p1<-ggplot(data=a, aes(x=NMDS1, y=NMDS2, color=smoking_current) )+
  geom_point()+
  stat_ellipse()+
  labs(y="MDS2", x="MDS1")+
  scale_shape_manual(values=c(16,17),
                     guide =
                         guide_legend(label.theme = element_text(size=6), 
                                      title=NULL)
                     )
p2<-ggplot(data=a, aes(x=NMDS1, y=NMDS3, color = smoking_current) )+
  geom_point()+
  stat_ellipse()+
  labs(y="MDS3", x="MDS1")+
  scale_shape_manual(values=c(16,17),
                     guide =
                         guide_legend(label.theme = element_text(size=6),
                                    title=NULL)
                     )

p1 + p2 + plot_layout(guides='collect') + plot_annotation(title="Current smoking") & theme_qk() & scale_color_uchicago()

Try some visualizations of pairs of factors using MDS1 & MDS2 (while 3D is preferred, stress of 3D is around .15, vs 0.22, so not huge difference) Separate HIV+/-, then use symbol colour to separate second factor.

p1<-ggplot(data=subset(a, candida=="Yes"), aes(x=NMDS1, y=NMDS2, color=HIV) )+
  geom_point()+
  stat_ellipse()+
  labs(y="MDS2", x="MDS1",
       title="Candida")+
  scale_shape_manual(
                     guide =
                         guide_legend(label.theme = element_text(size=6), 
                                      title=NULL)
                     )
p2<-ggplot(data=subset(a, candida =="No"), aes(x=NMDS1, y=NMDS2, color = HIV) )+
  geom_point()+
  stat_ellipse()+
  labs(y="MDS2", x="MDS1",
       title = "No Candida")+
  scale_shape_manual(
                     guide =
                         guide_legend(label.theme = element_text(size=6),
                                    title=NULL)
                     )
p_combined <- p1 + p2 + plot_layout(guides='collect') & theme_qk() & scale_color_uchicago()

p_ranges_x <- c(ggplot_build(p_combined[[1]])$layout$panel_scales_x[[1]]$range$range,
  ggplot_build(p_combined[[2]])$layout$panel_scales_x[[1]]$range$range)

p_ranges_y <- c(ggplot_build(p_combined[[1]])$layout$panel_scales_y[[1]]$range$range,
                ggplot_build(p_combined[[2]])$layout$panel_scales_y[[1]]$range$range)
p_combined & 
  xlim(min(p_ranges_x), max(p_ranges_x)) & 
  ylim(min(p_ranges_y), max(p_ranges_y))

p1<-ggplot(data=subset(a, HIV=="Yes"), aes(x=NMDS1, y=NMDS2, color=candida) )+
  geom_point()+
  stat_ellipse()+
  labs(y="MDS2", x="MDS1",
       title="HIV+")+
  scale_shape_manual(
                     guide =
                         guide_legend(label.theme = element_text(size=6), 
                                      title=NULL)
                     )
p2<-ggplot(data=subset(a, HIV =="No"), aes(x=NMDS1, y=NMDS2, color = candida) )+
  geom_point()+
  stat_ellipse()+
  labs(y="MDS2", x="MDS1",
       title = "HIV-")+
  scale_shape_manual(
                     guide =
                         guide_legend(label.theme = element_text(size=6),
                                    title=NULL)
                     )
p_combined <- p1 + p2 + plot_layout(guides='collect') & theme_qk() & scale_color_uchicago()

p_ranges_x <- c(ggplot_build(p_combined[[1]])$layout$panel_scales_x[[1]]$range$range,
  ggplot_build(p_combined[[2]])$layout$panel_scales_x[[1]]$range$range)

p_ranges_y <- c(ggplot_build(p_combined[[1]])$layout$panel_scales_y[[1]]$range$range,
                ggplot_build(p_combined[[2]])$layout$panel_scales_y[[1]]$range$range)
p_combined & 
  xlim(min(p_ranges_x), max(p_ranges_x)) & 
  ylim(min(p_ranges_y), max(p_ranges_y))

p1<-ggplot(data=subset(a, HIV=="Yes"), aes(x=NMDS1, y=NMDS2, color = smoking_current) )+
  geom_point()+
  stat_ellipse()+
  labs(y="MDS2", x="MDS1",
       title="HIV+")+
  scale_shape_manual(
                     guide =
                         guide_legend(label.theme = element_text(size=6), 
                                      title=NULL)
                     )
p2<-ggplot(data=subset(a, HIV =="No"), aes(x=NMDS1, y=NMDS2, color = smoking_current) )+
  geom_point()+
  stat_ellipse()+
  labs(y="MDS2", x="MDS1",
       title = "HIV-")+
  scale_shape_manual(
                     guide =
                         guide_legend(label.theme = element_text(size=6),
                                    title=NULL)
                     )
p_combined <- p1 + p2 + plot_layout(guides='collect') & theme_qk() & scale_color_uchicago()

p_ranges_x <- c(ggplot_build(p_combined[[1]])$layout$panel_scales_x[[1]]$range$range,
  ggplot_build(p_combined[[2]])$layout$panel_scales_x[[1]]$range$range)

p_ranges_y <- c(ggplot_build(p_combined[[1]])$layout$panel_scales_y[[1]]$range$range,
                ggplot_build(p_combined[[2]])$layout$panel_scales_y[[1]]$range$range)
p_combined & 
  xlim(min(p_ranges_x), max(p_ranges_x)) & 
  ylim(min(p_ranges_y), max(p_ranges_y))

References

Dixon, Kelly M., Geoffrey J. Cary, Graeme L. Worboys, and Philip Gibbons. 2018. “The Disproportionate Importance of Long‐unburned Forests and Woodlands for Reptiles.” Ecology and Evolution 8 (22): 10952–63. https://doi.org/gfs587.

Griffen, Ann L., Zachary A. Thompson, Clifford J. Beall, Elizabeth A. Lilly, Carolina Granada, Kelly D. Treas, Kenneth R. DuBois, et al. 2019. “Significant Effect of HIV/HAART on Oral Microbiota Using Multivariate Analysis.” Scientific Reports 9 (1): 19946. https://doi.org/gsfpz3.

Hutto, David, and Kyle Barrett. 2021. “Do Urban Open Spaces Provide Refugia for Frogs in Urban Environments?” Plos One 16 (1). https://doi.org/gr2r9n.